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ABSTRACT 



We investigate the fast (type III) migration regime of high-mass protoplanets or- 
biting in protoplanetary disks. This type of migration is dominated by corotational 
torques. We study the the details of flow structure in the planet's vicinity, the depen- 
dence of migration rate on the adopted disc model, and the numerical convergence of 
models (independence of certain numerical parameters such as gravitational soften- 
ing). 

We use two-dimensional hydrodynamical simulations with adaptive mesh refine- 
ment, based on the FLASH code with improved time-stepping scheme. We perform 
global disk simulations with sufficient resolution close to the planet, which is allowed 
to freely move throughout the grid. We employ a new type of equation of state in 
which the gas temperature depends on both the distance to the star and planet, and 
a simplified correction for self-gravity of the circumplanetary gas. 

We find that the migration rate in the type III migration regime depends strongly 
on the gas dynamics inside the Hill sphere (Roche lobe of the planet) which, in turn, 
is sensitive to the aspect ratio of the circumplanetary disc. Furthermore, corrections 
due to the gas self-gravity are necessary to reduce numerical artifacts that act against 
rapid planet migration. Reliable numerical studies of Type III migration thus require 
consideration of both the thermal and the self-gravity corrections, as well as a sufficient 
spatial resolution and the calculation of disk-planet attraction both inside and outside 
the Hill sphere. With this proviso, we find Type HI migration to be a robust mode 
of migration, astrophysically promising because of a speed much faster than in the 
previously studied modes of migration. 

Key words: accretion, accretion discs - hydrodynamics - methods: numerical - 
planets and satellites: formation 



1 INTRODUCTION 

In the standard model of planetary system formation giant 
planets form through core accretion outside the ice conden- 
sation boundary (4 — 5 AU), where the availability of water 
ice allows planetary cores to rapidly reach the critical mass 
of ^10 Earth masses, beyond which substantial gas accre- 
tion can occur ([Pollack et al.| [l996). This theory provides a 
good fit to the structure of the Solar System, but does not 
explain the presence of so-called 'hot Jupiters' (objects with 
minimum masses similar to or larger than Jupiter's mass, 
M7|_, and semi-major ax es a < 0.1 AU) discover ed in extra- 
solar planetary systeni s (|Mavor Quelo z 1995, Marcv et al.l 
l200Ql . lVogt et al.ll2002l ). Since the in situ formation of these 
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objects is difficult both in the core acc retion scen ario and 
through direct gravitational instability (Boss"2001), the in- 
ward migration of bodies and eccentricity pumping due to 
planet-planet and planet-remnant disc interaction become 
important parts of the new theory of planet formation. 

The standard t heory of p lanet miration 
("Goldreich k Tremaine "l979|, I198QI : | Lin Papaloizoul 
1993; Lin et al. 2000) considers the Lindblad resonances 
and neglects the corotational resonances, assuming a 
smooth initial density profile in the disc and a small density 
gradient at the corotation radius. Tidal torques due to the 
Lindblad resonances are believed to cause inward migration 
of a planet (type I or type II, for low and high mass planets 
respectively) from the region where it formed to the inner 
regions where many exoplanets are observed. The type II 
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migration times for giant planets are essentially viscous 
time-scales of disks, thus relatively long (|Wardlll997l V 



However, it was recently found that the co rotational res- 
onance can modify the type I migratio n mode (|Masset et al.l 
I2OO6I . IPaardekooper Mellemalbooih or lead to a new and 
very fast migration mode (type III, or run- away migra- 
tion) that depends strongly on the gas flow in the planet's 
vicinity and does not have a predetermined direction 
(iMasset Papaloizoul l2003l : lArtvmowic'3 l2004l : IPapaloizoul 
I2OO 5': 'Artvmowicz"2QQ6') . This type of migration was studied 
numerically by Masset & Papaloizou (2003) who performed 
two-dimensional simulations of a freely migrating planet and 
a steady state migration with fixed migration rate d for a 
range of the migration rates. Global, high resolution two- 
and three-dimensional simulations of freely migrating plan- 
e ts w ere performed by D'Angelo et al. (2005). Papaloizou. 
(j2005l l considered local shearing box simulations. 



These numerical simulations showed that the planet's 
orbital evolution depends strongly on th e choice of the 
simulation parameters, e.g. grid resolution (D'Angel o et al.l 
I2OO5I ). softening of the planet gravitation potential, etc. The 
reason for this dependence can be the simplifications com- 
monly used in the disc model. The most important simpli- 
fications are the use of two-dimensional simulations, ignor- 
ing self-gravity and the local-isothermal approximation, that 
imposes a static temperature distribution in the planet's 
vicinity. 

This is the first in a series of papers devoted to type 
III migration of high-mass protoplanets interacting with the 
protoplanetary disc. In this paper we study the dependence 
of the planet migration on the applied disc model and we 
defer the description of the ph ysics of type III migration it - 
self to Peplinski et al] (|2007al ) and IPephnski etM (|2007bl ) 
(henceforth Paper II and Paper III). We concentrate on two 
aspects: a modification of the local-isothermal approxima- 
tion, that allows an increase the temperature inside the 
Roche lobe, and a correction of the gas acceleration (due 
to the gas self-gravity) that forces the circumplanetary disc 
to move together with the planet. 

In order to analyse this we study the gas flow in the 
planet's vicinity by performing numerical hydrodynamical 
simulations in two dimensions of a gaseous disc interacting 
with the star and one planet. The planet is allowed to mi- 
grate due to disc-planet interaction and we can study its 
orbital evolution in the non-steady state. Since type III mi- 
gration is so sensitive to the gas flow near the planet, we 
include a careful analysis of the effects of various numerical 
parameters which influence these flow patterns. 

The layout of the paper is as follows. In Sections [2] and 
[3] we give the basic equations, describe the disc model and 
the numerical method. Sections |4] and [5] contain the conver- 
gence tests and a discussion of the effects of various modifi- 
cations of disc model and numerical parameters for inward 
and outward migration. Finally in Section [6] we discuss their 
implications for numerical simulations of type III migration. 



2 DESCRIPTION OF THE PHYSICAL MODEL 
2.1 Disc model 

We adopt in our simulations a two-dimensional, infinitesi- 
mally thin disc model and use vertically averaged quantities, 
such as the surface mass density 



00 

I 



pdz, 



(1) 



where p is the mass density. We work in the inertial refer- 
ence frame, in a Cartesian coordinate system (x,y,^), also 
see Sect. [321 The plane of the disc and the star-planet sys- 
tem coincides with the z — plane. The centre of mass of 
the star- planet system is initially set at the origin of the co- 
ordinate system. Since the star and the planet are allowed to 
migrate freely due to the gravitational interaction with the 
disc, and the total angular momentum of the whole system 
is not fully conserved due to the open boundary conditions 
(see Section [2. 3. 2 1) . the centre of mass can move slowly away 
from the origin of the coordinate system. The positions and 
masses of the star and the planet we denote by rs, Ms, rp 
and Mp respectively. 

The gas in the disc is taken to be inviscid and non-self- 
gravitating. The evolution of the disc is given by the two- 
dimensional {x,y) continuity equation for S and the Euler 
equations for the velocity components v = (vx^Vy). These 
equations can be written in conservative form as 

-hV(Si;) = 0, (2) 



dt 



dt 



(3) 



where P is two-dimensional (vertically integrated) pressure, 
and $ is the gravitational potential generated by protostar 
(subscript S) and planet (subscript P) 



GMs 
r -rs\ 
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■ rp 



(4) 



We do not consider the energy equation, since we use a 
local isothermal approximation (see Sect. [2XT|) . 

The gravitational potential close to the star and the 
planet is softened in the following way: 



-(1.875^^ -7^^ + 7.875^^+ 
-4.375^2 + 2.625) 



for ^ > 1 



for ^ ^ 1 



(5) 



where ^ = r/rsoft and rsoft is the so-called smoothing length 
(or gravitational softening). The reason for using this for- 
mula instead of the standard one (where r and rsoft are 
added geometrically) is to remove the dependence of the Ke- 
plerian speed on the gravity smoothing length for r > rsoft 
in the disc surrounding the body. This is especially impor- 
tant for simulations using a Cartesian grid since the stel- 
lar gravity has to be softened too. Moreover it allows for 
a bigger smoothing length around the planet without influ- 
encing the outer part of the Roche lobe. This is necessary 
because the planet is moving across the grid and its position 
with respect to cell centres varies. As Nelson Beiiz (200^al ) 
pointed out, in this case, in order to avoid unphysical effects 
on the planet's trajectory caused by close encounters with 
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Figure 1. The real disc aspect ratio calculated with respect to 
the planet for two formulae for isothermal sound speed Cg. Curve 
1: EOS2 (Eq.[n]); curve 2: EOSl (Eq.[lO]). Curve 3 shows the real 
aspect ratio for EOS2 after taking into account the gravitational 
softening (the result for EOSl is not affected by this). 



cell centres, the smoothing length needs to be larger than 
half the dimension of the cell. In our simulations rsoft is at 
least a few times the cell size. 

Unlike the standard formula for softening of the gravity, 
our formula corresponds to a spherical body with a finite 
radius given by rsoft- The internal structure of this body 
can be found using the Poisson equation and is given by the 
formula 



; for ^ > 1 



^soft 



-(-78.75^^ + 210^3+ 

-157.5^^+ 26.25) ; for ^ ^ 1 



(6) 



Since r = rsoft is the "surface" of the body generating 
the gravitational potential, we will include the mass of the 
gas within the gravitational softening Msoft to the effective 
planet mass. For more discussion see Sections 12.1.21 , 15.1.21 
and [523] 

Previous studies of migrating Jupiter-mass planets 
showed that an important issue is the tr eatment of the 
torques arising from within the Hill sphere ([Nelson &: Bend 
l2003bl : lD'Angelo et al.ll2005l ). It is often assumed that these 
torques are strong, but nearly cancel, and thus can be ne- 
glected. In addition it has been argued that accelerating a 
planet through its own (bound) envelope would be unphys- 
ical. In general, the exclusion of the torques arising from 
within the Hill sphere is inconsistent, but can be justified 
if the density distribution in the Hill sphere is (almost) 
static and if the circumplanetary disc can be considered as 
a separate system. As we show below, these assumptions 
are not necessarily satisfied in the highly dynamic case of 
type HI migration. The torques from the planet's vicinity 
need be taken into account. Obviously, these torques de- 
pend strongly on the density distribution near the planet, 
which in a numerical simulation will depend on a combina- 
tion of the resolution and the parameters determining this 
density distribution. The gravitational softening rsoft is one 
of these parameters. In Sect. 15.1.11 we study the effects of 
gravitational softening. 



2.1.1 Equation of state 

Since a self-consistent calculation of the temperature is pro- 
hibitively expensive, we adopt the usual local isothermal ap- 
proximation: the disc is treated as a system having a fixed 
temperature distribution. In this case the equation of state 
has the form 



(7) 



Cs being the local isothermal sound speed. The thermal state 
of the fluid is given by Cs, which is usually assumed to be a 
power law of the distance to the star rs = |r — rs|. Assuming 
that the material inside a circumstellar disc is in hydrostatic 
equilibrium and neglecting the planet gravitational field, we 
can get the simple, widely used formula 



Cs - HQs, 

where H is the disc scale height and 
angular velocity in the circumstellar disc 



Qs 



GMs 



(8) 

is the Keplerian 
(9) 



This formula, for the constant disc opening angle, gives 

Cs = hsTs^s, (10) 

where hs = H/rs is the disc aspect ratio with respect to the 
star. In this case the Mach number of the flow in the disc 
is a constant given by 1 //is • We will designate this choice as 
EOSl (equation of state 1). 

For EOSl the gas temperature is not modified by the 
presence of the planet. Because of the gravitational field of 
the embedded planet, the disc aspect ratio calculated with 
respect to the planet hp = H/\r — rp| can achieve very 
small values (see Fig.[T]), allowing the planet to collect large 
amounts of material. Physically, we should expect a tem- 
perature increase in the vici nity of a Jupiter - mass planet. 
Two-dimensional simulations (|D'Angelo et al.i r2Q03) showed 
that an accreting Jupiter on a constant orbit can have a cir- 
cumplanetary disc with constant hp ranging from 0.2 to 0.4, 
and three-dimension al simulations show even higher values 
felah r k Ki^l2QQ6h . In the highly dynamic case of a rapidly 
migrating planet hp can be expected to attain similar val- 
ues. For this reason we introduce a prescription for the sound 
speed, which depends on the distance to both the star and 
the planet 



((/isrs)^ + {hprp)^y/ 



(11) 



where n is a non-dimensional parameter, rp = |r — rp| is 
the distance to the planet, and Qp is the Keplerian angular 
velocity in the circumplanetary disc 



GMp 



(12) 



This equation gives a constant disc aspect ratio in the cir- 
cumstellar disc hs far away from the planet, and a con- 
stant disc aspect ratio in the circumplanetary disc hp in the 
planet's vicinity. Parameter n = 3.5 is chosen to smoothly 
join equations ()11|) and (|lQp . We will refer to this approach 
as E0S2. 

As pointed out above, the torques from within the Hill 
sphere are important for type HI migration, and will depend 
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on the choice for hp (as well as rsoft)- Tests of how hp in- 
fluences the convergence behaviour of our simulations are 
presented in Section l4?2l 



2.1.2 Corrections for the gas self-gravity 

For low mass discs, it is usual to not include the effects of 
gas self- gravity in numerical simulations. The argument is 
that these effects are minor in not too massive disks, while 
the calculation of self- gravity is particularly expensive. How- 
ever, in the case of type III migration the planet can collect 
a considerable envelope with a mass equal to its own mass. 
The planet migrates through interaction with the disc ma- 
terial, but without self-gravity, its envelope will not do the 
same. This may lead to an artificial increase of the planet's 
inertia, due to the discrepancy between the planet position 
and the centre of mass of the planet's gaseous envelope, and 
can result in an a dditional, non-physical f orce acting against 
pla net migration (jPapaloizou et all l2007l l . 

iD'Angelo et al.1 moi )' investigated type III migration 
and found that increasing the spatial resolution of their 
simulations led to gradually slower migration. Since higher 
resolution allows for more mass to accumulate within the 
planet's Hill sphere, we believe that their results reflect ex- 
actly the effect described above. 

Instead of calculating the full effect of self-gravity, we 
apply a first order correction by forcing the planet and its 
gaseous envelope to move together. This is done by modi- 
fying the acceleration of the gas in a planet's vicinity (ag), 
adding the planet acceleration due to the torque from the 
gas ax to the acceleration exerted on the gas by the star 
and the planet: 



GMs{r - rs) GMp{r - r?) 



+ 



+aT max(0, 1 



(rp/re: 



(13) 



where rs, tp are the positions of the star and the planet, 
respectively. Mp is taken to be either the planet mass Mp, 
or Mp, the planet mass plus all the gas mass within the 
smoothing length rsoft around the planet. The latter is equiv- 
alent to assuming that the gas inside rsoft attains the spher- 
ical density distribution given by Eq. (|6]). 

The last term gives the planet acceleration multiplied 
by a function limiting the correction to the planet's enve- 
lope using a parameter renv renv should be of order the 
size of the circumplanetary disc with Ru > renv > rsoft- 
We find that renv — 0.5-Rh removes the artificial increase of 
the planet's inertia, provided the density distribution inside 
renv is relatively symmetric and smooth, as is the case in 
our simulations. 

Notice that renv and rsoft are independent parameters, 
^soft gives the position of the "surface" of the protoplanet 
(radius where $ starts to differ from the point-mass poten- 
tial), whereas renv gives the size of the region that dynam- 
ically belongs to the planet and should follow the planet 
in its radial motion. Our correction for the gas accelera- 
tion allows to reduce the non-physical eccentricity of the 
orbits in the circumplanetary disc driven by the planet's ra- 
dial motion. The renv is thus a third parameter influencing 
the density distribution near the planet (the others being 
the gravitational smoothing, rsoft, and the disc aspect ra- 



tio in the circumplanetary disk, hp). The effects of renv on 
the numerical convergence for type HI migration are de- 
scribed in sections 14.11 and 15^2.11 The differences between 
using Mp = Mp and Mp = Mp are studied in Sect. 15.1.21 

We stress that the described method is a crude approx- 
imation introduced to remove the non-physical effects from 
the planet's orbital evolution. It does not allow any detailed 
study of the real flow of self- gravitating gas inside the Roche 
lobe or of the structure of the planet's gaseous envelope, 
which are beyond the scope of this paper. 



2.1.3 Accretion onto the planet 

The orbital evolution can also be affected by gas accretion 
onto the planet. This can be dealt with either by using Mp — 
Mp instead of Mp in Eq. [131 or by removing matter from 
the planet's neighbourhood |r — rp| < race, and adding it 
to Mp. In the latter case, we also add the momentum of the 
accreted gas to that of the planet's. When removing gas from 
the disc, this is done after each integration step according 
to the formula 



■ rp 



(S - S) , (14) 



AS — max ( 0, max ( 0, 1 

"Tacc 



where race is the accretion time-scale and S is an average sur- 
face density in the region defined by race < Ir — rp| < 2racc. 
The size of the accretion region race is defined by the size of 
the smoothing length rsoft, i-e. the size of the gravitational 
source. In simulations where we explicitly include accretion, 
we use 



0.5rs. 



(15) 



which is an order of magnitude smaller than the Roche lobe 
size. The effects of the different ways to deal with accretion 
are studied in Sect. 15.1.21 



2.2 Equation of motion for the star and the planet 

The goal of the current study is to investigate the orbital 
evolution of the planetary system due to the gravitational 
action of the disc material. Since the calculations are done 
in the inert ial reference frame, the equation of motion for 
the star and the planet have the simple form: 



rs 



GMp{rs-rp) 
\rp - rsf 

GMsjrp - rs) 
Irp - rsf 



■I 
■/ 



G{rs -r)dMB{r) 
\r - rsf 

G{rp - r) dMuir) 



■ rp 



(16) 



(17) 



In both cases the integration is carried out over the disc 
mass Md included inside the radius raise (thus removing 
the corners of the Cartesian grid from consideration). 



2.3 Simulation setup 

In the simulations we adopt non-dimensional units, where 
the sum of star and planet mass Ms + Mp represents the 
unit of mass. The time unit and the length unit are chosen 
to make the gravitational constant C = 1. This makes the 
orbital period of Keplerian rotation at a radius a = 1 around 



© 0000 RAS, MNRAS OOO.fUfTT] 



Type III migration; Disc model 5 



\2 
\ \. 3 



4 '\ 



0.0035 
0.003 
0.0025 
0.002 
0.0015 
0.001 
0.0005 


0.5 1 1.5 2 2.5 3 3.5 4 
r 

Figure 2. The initial surface density profile for oly. ranging from 
0.0 to —1.5 (outward migration case). 



a unit mass body equal to 27r. However, when it is necessary 
to convert quantities into physical units, we use a Solar-mass 
protostar Ms = ^ Jupiter-mass protoplanet Mp = 

Mt^, and a length unit of 5.2 AU. This makes the time unit 
equal to 11.8/27r years. 

In all the simulations the grid extends from —4.0 to 4.0 
in both directions around the star and planet mass centre. 
This corresponds to a disc region with a physical radius of 
20.8 AU. 



The smoothing length of the stellar potential, rsoft, is 
taken to be 0.5. Unless otherwise noted, for the planet the 
standard value is 0.33i?H, where Ryl = a[Mp/(3Ms)]^^/^\ 
the Hill radius. The corresponding size of the envelope renv 
in Eq. ([13]) is set to rsoft or 0.5i?H. 

If accretion onto the planet is included, we use an ac- 
cretion time-scale race = IOtt and race = 0.5rsoft- 



2.S.2 Boundary conditions 

To stop reflections off the boundary, we use an outflow- inflow 
boundary condition with a so-called killing wave zone. This 
zone extends from radius 3.65 to 3.9. In this region the solu- 
tion of the Euler equations for X {X stands for the surface 
density or the velocity) is changed after each time-step by 

AX=-{X-Xo)ct>{r), (20) 

where Xq is the initial condition (disc in sub-Keplerian ro- 
tation), Td is the damping time-scale, and (j){r) is a function 
increasing from at the inner boundary of the killing wave 
zone to 1 at the outer boundary of the killing wave zone. 
Outside the radius 3.9 the solution X is replaced by Xq. In 
our simulations we use Td = 18. 

On our Cartesian grid the disc cannot be treated as an 
isolated system and some mass and angular momentum flow 
through the boundary is unavoidable. However, the losses 
are usually relatively small. 



2.3.1 Initial conditions 

The initial surface density S profile is given by a modified 
power law: 



: ^(rc)i;o(rc/ro)'' 



(18) 



where rc = Ir — rc| is the distance to the mass centre of the 
planet-star system, ro is a unit distance, and '0 is a function 
that allows introducing a sharp edges in the disc (see Fig. [2]) • 
We characterize the disc mass by the disc to the primary 
mass ratio 

In the simulations /xd ranges from 0.001 to 0.01. For the 
Minimum Mass Solar Nebula (MMSN) /xd = 0.00144 (for 
as = —3/2). We investigate different density profiles by 
changing from —1.5 to 0.0. 

Since we focus on type HI migration only and do not 
analyse the problem of orbital stability inside a gap, we do 
not introduce the planet smoothly nor keep it on a constant 
orbit for the time needed to create a gap. For most cases 
with as < 0.0 the density gradient given by the initial pro- 
file is sufficient to start rapid inward migration. In the case 
of outward migration we start migration by introducing an 
additional density jump at the planet's position. 

The planet mass is taken to be Mp/Ms = 0.001 (i.e. one 
Jupiter mass, Mt^, for a one-solar-mass star). The planet 
starts on a circular orbit of semi-major axis equal 3.0 and 
0.8 for the inward and outward migration case respectively. 

The aspect ratio for the disc with respect to the star is 
fixed at hs = 0.05, whereas the circumplanetary disc aspect 
ratio hp is taken from the range 0.2 to 0.6. 



3 DESCRIPTION OF THE NUMERICAL 
METHOD 

3.1 Code 

We adopted the FLASH hydro code version 2.3 written by 
the FLASH Code Group from the Centre for Astrophysical 
Thermonuclear Flashes at the University of Chicago in 
1997. 

FLASH is a modular, adaptive- mesh, parallel simula- 
tion code capable of handling general compressible flow 
problems. It is designed to allow users to configure ini- 
tial and boundary conditions, change algorithms, and add 
new physics modules . It uses the PARAMESH library 
(jMacNeice et al.l l200Q') to manage a block-structured adap- 
tive mesh, placing resolution elements only where they are 
needed most. PARAMESH consists of a set of subroutines 
which handle refinement /derefinement, distribution of work 
to processors, guard cell filling, and flux conservation. It also 
uses the Message- Passing Interface (MPI) library to achieve 
portability and scalability on parallel computers. For our 
purpose, the code is used in the pure hydrodynamical mode 
in two dimensions, and the adaptive mesh is used to achieve 
high resolution around the planet (see Sect. 13. 2|) . 

Euler's equations are solved using a directionally 
split version of the Piec ewise- Parabolic Method {PPM, 
IColella &: Woodward! 1 1983 ). It represents the flow variables 
inside the cell with piecewise parabolic functions. This 
method is considerably more accurate and efficient than 
most formally second-order algorithms. 



http:/ /flash. uchicago.edu 
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Figure 3. Convergence test for the inward migrating Jupiter. The left panel shows the time evolution of the planet's semi-major axis a, 
the right panel the mass of the gas inside a Hill sphere • Curves 1 and 2 correspond to simulations with 4 and 5 levels of refinement 
and Tenv = (planet acceleration is not added to the gas acceleration); curves 3 and 4 correspond to simulations with 4 and 5 level of 
refinement and renv = 0.5i?H- The time unit is the orbital period of the body circulating at a = 1. The planet mass Mp = 0.001. 



To solve the equations of motion for the star and the 
planet we use a 4-th order Runge-Kutta method. To keep 
the second order accuracy in time PPM needs information 
about the gravitational field from the beginning and the end 
of the time- step. Therefore we perform the integration of the 
planet's orbit in two steps. First we integrate using the force 
exerted by gas taken from the previous hydro- step, This 
provides PPM with the approximate position of the planet 
at the end of the current time- step. Second, after updating 
the hydrodynamic quantities, we correct the planet position 
using a linear interpolation between the old and new values 
of acceleration. 



3.2 IMesh and grid structure 

High resolution is needed to calculate accurately the gas flow 
in the planet's vicinity. For this purpose we use the Adap- 
tive Mesh Refinement (AMR) module included in FLASH. 
AMR allows us to change the grid structure during the sim- 
ulation, and gives the possibility for the refined region to 
follow the planet's motion. The resolution is increased by a 
factor of 2 between refinement levels. Our simulations use 
a lowest resolution mesh of 800 cells in each direction, and 
a square region around the planet is refined. The maximal 
cell size in the disc (lowest level of refinement) is about 1% 
of the smallest value of the planet's semi-major axis. The 
cell size close to planet is at least 1.8% of Hill sphere radius 
(corresponding to 4 levels of refinement). 

We use a Cartesian grid for our calculations. This choice 
was made because the more usual cylindrical (co-rotating) 
grid geometry has few benefits for the case of a rapidly 
migrating planet, and also suffers from variable (physical) 
resolution with radius. However, this does not mean that a 
Cartesian grid is problem free. The most important problem 
is the diffusivity of the code visible in Ide Val-Borro et aP 
(l2006lV however the time-scale of this process is relatively 
long and does not influence the rapid migration. It is dis- 
cussed in Appendix lAl 



3.3 IVLulti-level time integration 

In order to reduce the computational time we extended 
FLASH with the option of different time- steps for differ- 
ent refinement levels. The algorith m used is very similar t o 
nested-grid technique presented in iD'Angelo et all (|2QQ2l ). 
Coarser blocks can use longer time- steps than finer blocks, 
but all blocks at a given resolution have to use the same 
time- step. In this case a time- step is a monotonic function 
of the grid resolution (in FLASH all time- steps have to be 
power of 2 multiples of the time-step being used on the finest 
block). This resolution dependent time-step can save a lot 
of computer time, but perhaps more importantly reduces 
the numerical diffusion in blocks at coarser refinement lev- 
els, which need a fewer number of time-steps to calculate a 
given evolution time. 



4 NUIVLERICAL CONVERGENCE 

The problem of type HI migration is numerically challenging 
since we have to deal with the complex evolution of the 
co-orbital flow in a partially gap opening regime, close to 
the planet. Clearly one can expect the results to depend on 
the disk model close to the planet, and in this section we 
investigate how the convergence behaviour (dependence on 
resolution) depends on our choices for the correction for self- 
gravity (renv), circumplanetary disc aspect ratio (hp), and 
gravitational softening (rsoft)- Here we discuss the inward 
migration case only, but the conclusion are valid for the 
out ward directed migratio n too. 

iD'Angelo et al.l (|2005h showed that it is hard to achieve 
numerical convergence for these types of flows. They found 
the migration to be highly dependent on the torque calcu- 
lation prescription, and on the mesh resolution (see their 
Fig. 8 and 9). We can qualitatively reproduce their results 
when using EOSl, and for renv = and Mp = Mp (Eq.[l3|), 
see curves 1 and 2 in Fig. [S] In this case higher resolution 
allows more mass to accumulate in the planet's vicinity, and 
the lack of any self-gravity allows this mass to increase the 
planet's inertia, slowing down its migration. 
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Figure 4. Surface-density radial cut through the Roche lobe at 
the planet position for different temperature profiles. Curve 1 cor- 
responds to EOSl (Eq. JlO]); curves 2, 3 and 4 to EOS2 (Eq. ^ 
with hp equal 0.2, 0.3 and 0.4 respectively. 

4.1 Correction for gas self-gravity 

As explained in Sect . 12 . lT2l neglecting self- gravity is likely to 
be one cause of this bad convergence behaviour. To explore 
this further we use the correction for self-gravity described 
in Sect. 12.1.21 The key parameter of this correction 
the size of the region where the gas is dynamically connected 
to the planet. To establish its value, we performed a series of 
simulations with different resolutions (3, 4, and 5 refinement 
levels) and different values of renv' 0.0 (no correction), rsoft 
(the smoothing length for gravity) and 0.5i?H- The last value 
is bigger than the estimated radius of the circumplanetary 
disc. From the simulations we know that the circumplan- 
etary disc extends up (0.3 — 0.4)i?H (see Fig. [5])- We can 
expect Tenv to bc bigger than the radius of the circumplan- 
etary disc, since Eq. ([13]) does not define a sharp edge to 
the planet's envelope. Some of the results of these tests are 
presented in Fig. [3] In these simulations we used a constant 
Mp Mp. 

We see that the correction reduces the 'additional iner- 
tia' and increases the migration rate, even though the mass 
of the gas accumulated in the Hill sphere increases signifi- 
cantly (right panel). For renv = Tsoft the correction region is 
smaller than the real size of the planet envelope, and we find 
results similar to curves 1 and 2. The convergence is much 
better for renv = O.bRu (curves 3 and 4), but even there we 
have not reached complete convergence with 5 refinement 
levels. This is caused by the mass accumulation in the cir- 
cumplanetary disc. Comparing the simulations with 4 and 
5 refinement levels (curves 3 and 4 respectively) we can see 
that the planet migrates faster after the mass accumulation 
in the circumplanetary disc stops. As the amount of mass 
accumulation depends on the assumed thermal structure of 
the circumplanetary disc (here taken to follow EOSl), we 
will explore this point more in the next section. 

One may consider excluding the Roche lobe interior 
from the torque calc ulation as another solution t o the artifi- 
cial inertia problem (Masset Pa paloizoul 12003). However, 
the gas flow in the planet vicinity is very variable and it is 
impossible to define a single cutting radius to remove the 
region dynamically connected to the planet from the torque 
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Figure 5. The flow in the planet's proximity for the rapid inward 
migration. In this case the flow symmetry is broken and the big 
ammount of the mass is accumulated in the planet's proximity. 
The plot shows the surface density distribution and the flow lines 
in the frame co-moving with the planet and covers a square region 
of the size of 4i?H- The colour scale is logarithmic. The pink lines 
show the approximate borders of the different regions in the disc: 
inner disc I, outer disc II, horseshoe region III, circumplanetary 
disc IV, co-orbital flow transferring the gas from the inner to the 
outer disc V and the gas stream entering the circumplanetary disc 
VI. The arrows show the direction of the flow. 

calculation and at the same time keep all the flow-lines that 
belong to the corotational flow. Moreover, in the case of the 
fast migration, the interior of the Roche lobe is not sepa- 
rated from the corotational flow and we cannot neglect the 
gas inflow into the circumplanetary disc (Fig. [5])- 

4.2 Temperature profile in the circumplanetary 
disc 

As explained in Sect. 12.17X1 the use of EOSl imposes a very 
thin circumplanetary disc. In such a thin disc waves can 
freely propagate and the shocks in the circumplanetary disc 
can end very close to the planet (see Fig. 2] curve 1). Such 
a strong shocks modify the flow inside the Roche lobe, redi- 
recting the gas to the planet's proximity. In the case of a 
planet on a constant orbit, this configuration does not cause 
any serious problems, since the co-orbital flow is weak and 
the gas inflow into the circumplanetary disc is limited. This 
changes for a migrating planet. In this case the inflow is 
potentially much larger since the planet moves into previ- 
ously undisturbed disk regions. In fact, mass is often seen 
accumulating in the circumplanetary disc. 

The flow structure of the gas in the planet's proximity 
for the rapid inward migrating planet is presented in Fig. \5\ 
On the plot we indicate the four important regions: the in- 
ner disc I, the outer disc II, co-orbital region HI and the 
circumplanetary disc IV. For slow (type II) migration the 
structure of the shocks in the circumplanetary disc and the 
flow pattern in the Hill sphere is almost symmetric, but for 
the rapidly migrating planet the flow symmetry is broken. 
The spiral shocks in the circumplanetary disc are very time 
dependent and are frequently replaced by a one-arm spiral. 
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Figure 6. Convergence test for the inward migrating Jupiter with the disk aspect ratio hp = 0.4. The left panel shows the time evolution 
of the planet's semi-major axis a, the right panel the mass of the gas inside a Hill sphere Mh- Curves 1, 2 and 3 correspond to 3, 4 and 
5 levels of refinement. The planet mass Mp = 0.001. 
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Figure 7. Surface density and the flow lines in the planet's vicinity for an inward migrating Jupiter. Different plots correspond to 
different temperature proflles. Upper-left: EOSl (no dependence on the planet's position). Upper-right, lower-left and lower-right: EOS2 
with hp is equal 0.2 , 0.3 and 0.4 respectively. The plotted domain is a square region of the size of 4i?H- The colour scale is logarithmic. 
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Figure 8. The planet's orbital evolution for four discussed models: the simplest model without any corrections (curve 1), the model 
with correction of the gas acceleration (curve 2), the model with changed temperature profile (curve 3) and the final model with with 
the correction to the gas acceleration applied and the temperature rising in the planet vicinity (curve 4). The left panel shows the time 
evolution of the planet's semi-major axis a, the right panel the mass of the gas inside a Hill sphere Mh- The planet mass is Mp = 0.001. 



It is caused by a strong mass inflow into the circumplane- 
tary disc. The horseshoe region shrinks to a single tadpole- 
like region and a strong co-orbital flow (marked V and VI) 
develops. The co-orbital flow interacts with the bow shocks 
and finally is divided into a flow transferring matter from 
the inner disc into the outer disc and a stream entering the 
Roche lobe and accumulating in the circumplanetary disc. 

Due to the strong co-orbital flow a significant amount 
of material is moved into the Roche lobe and starts inter- 
acting with the shocks originating in the circumplanetary 
disc. Since the disc aspect ratio h decreases with decreasing 
distance to the planet, the shocks become stronger and the 
material moves along the shock directly to the central den- 
sity spike (see curve 1 on Fig.|4}. The radial size of this spike 
is given by gravitational softening rsoft, since at this posi- 
tion h grows and the shocks disappear. Inside this region 
(about 0.6rsoft) gas moves on circular orbits. When using 
EOSl any amount of material can be added to this disc, 
since we neglect self-gravity and the gas does not heat up 
during compression. As a result the mass of the circumplan- 
etary disc is only limited by the grid resolution and the value 
of rsoft. 

This is the cause of the failing numerical convergence 
for the runs with renv = 0.5i?H described in the previous 
section. From Fig. [3] we see that the migration rate for the 
simulation with 4 refinement levels (curve 3) increases af- 
ter the inflow into the circumplanetary disk stops at about 
40 orbits. This happens because gas moved from the inner 
to the outer disc exchanges with the planet twice as much 
angular momentum as the 'accreted' gas. 

To address this problem we replace EOSl with E0S2. 
The latter maintains a constant disc aspect ratio in the cir- 
cumplanetary disc and keeps h independent of the gravita- 
tional softening. To test the effects of this different equation 
of state, we performed simulations with different aspect ra- 
tios of the circumplanetary disc (/ip equal 0.2, 0.3 and 0.4) 
at different resolutions (3, 4 and 5 refinement levels), while 
neglecting our self-gravity correction (renv = 0) and keeping 
the planet mass constant at Mp during the whole simula- 
tion. 



We find that the model with /ip — 0.2 gives no im- 
provement over the EOSl results, the case with /ip = 0.3 
gives some, but only for hp — 0.4 numerical convergence 
for 5 refinement levels was achieved. These results can be 
understood by studying the density profile near the planet 
(Fig. 131) and density and flow patterns (Fig.E]). The /ip = 0.2 
case gives a profile very similar to EOSl. The model with 
/ip = 0.3 results in a somewhat smoother profile, but only 
for hp — 0.4 the density profile is smooth enough for the 
shocks to become unimportant, and the gas inflow into the 
circumplanetary disc is stopped, resulting in a lower mass 
inside the Roche lobe. The planet's orbital evolution and the 
mass of the Hill sphere content for the /ip = 0.4 simulation 
at different refinement levels are shown in Fig. [6] 

These results seem to indicate that the change of tem- 
perature profile alone is sufficient for achieving numerical 
convergence, even if we do not apply any corrections for self- 
gravity. The reason for this is the smoother density distri- 
bution and the lower mass of the gas in the circumplanetary 
disc. The 'artificial inertia' effects are particularly strong 
when the mass of the circumplanetary disk is high, and its 
density profile strongly peaked. In this case any small dis- 
crepancy between the disc centre and planet position results 
in strong, high frequency oscillations of all orbit parame- 
tersB 

However, even for the changed temperature profile, 
some 'artificial inertia' effects remain. To illustrate this we 
compare the planet's orbital evolution for four different mod- 
els in Fig.[8l Model 1 is the simplest model (without any self- 
gravity correction and with EOSl), model 2 uses the self- 
gravity correction (renv = O.Si^n), model 3 employs E0S2 
(/ip = 0.4) but no self-gravity correction, and model 4 ap- 
plies both (renv = 0.5i?H and hp — 0.4). We see that models 
2, 3 and 4 show a higher migration rate than model 1. For 
model 2 there is a visible kink at 45 orbits, which is the mo- 
ment when the limiting mass of the circumplanetary disc is 
reached. This kink is invisible for model 4, since the 'accreted 

^ This oscillations are not visible on the plots, since we removed 
them by averaging over 5 orbits. 
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Figure 10. The torque exerted by the gas on the planet for simu- 
lations with different smoothing length of the planet gravitational 
field Tsoft equal 0.0208 (curves 1, 2 and 3) and O.SRu (curves 4, 
5 and 6). The torque was calculated separately for the interior of 
the planet's gravitational softening rp < rgoft (curves 1 and 4) 
and the rest of the Roche sphere rgoft < < Rh (curves 2 and 
5), where rp is the d istance to the planet. The curves 3 and 6 
give the torque from the whole Hill sphere Frl (the sum of the 
curves 1, 2 and 4, 5 respectively). 



mass' is an order of magnitude smaller, and this mass is al- 
ready reached after the few first orbits. Comparing models 3 
and 4 we see that the latter has a larger migration rate than 
the third model, implying that the change of the tempera- 
ture profile does not remove the effects of 'artificial inertia' 
completely. We hence should consider both mechanisms to 
be equally important. 



5 TORQUE CALCULATION AND EFFECTIVE 
PLANET MASS 

In the previous section we discussed the modifications of the 
disc model that allow to achieve the numerical convergence. 
In this chapter we will concentrate on the calculation of the 
torque inside the Hill sphere, and the models with the vary- 
ing planet mass. Especially we will focus on the dependence 
of the migration on the choice of Mp. We will discuss it for 
the inward and outward migration case separately. 



5.1 Inward migration case 

5.1.1 Gravitational softening 

When including all material inside the Roche lobe in the 
calculation of the torque, the value for the gravitational 
softening rsoft can play an important role, even though in- 
creasing the value of the circumplanetary disc aspect ratio 
hp should make a less dependent on rsoft- Indeed we find 
that for hp ^ 0.4 the migration behaviour does not depend 
on rsoft- The results of the runs with hp — 0.4 and rsoft 
equal 0.0208 and O.SRu are presented in Fig. [9] In the first 
case rsoft is constant during the simulation and its value in 
Rh units is ranging from O.lRu (initial value) up to about 
0.25i^H- In both simulations we used Mp = Mp. The first 
and the second plot show the planet orbital evolution and 



the mass of the gas inside a Hill sphere. The orbital evolu- 
tion is independent on the size of smoothing lengths of the 
planet's gravitational field, however a larger rsoft allows the 
planet to accumulate larger amounts of material. 

Above we have argued against excluding any part of 
the disc in the torque calculation. The results in this section 
allow this to be illustrated better. Figure [10] presents the 
torque exerted on the planet by the gas contained within 
the Hill sphere. Curves 3 and 6 show the torque from the 
entire Hill sphere Frl- We divided this torque into two parts, 
namely the contribution from rp < rsoft (Tsoft, curves 1 and 

4) and the contribution from rp > rsoft (Tout, curves 2 and 

5) , where rp is the distance to the planet. The two set of 
curves (1, 2, 3) and (4, 5, 6) correspond to rsoft equal 0.0208 
and O.SRu respectively. 

Frl is driving the migration during the fast migration 
phase (first 60 orbits) and drops in the slow, type II like 
migration phase. As we saw before, it is almost indepen- 
dent on rsoft- However, the partial torques Fsoft and Fout 
are varying strongly between the two choices for rsoft- Dur- 
ing the fast migration phase even the sign of Fout depends on 
rsoft- During the slow phase the two choices for the soften- 
ing/cutting radius both give a positive Fout, but its value in- 
creases with decreasing rsoft , consistent with what was found 
by IP'Angelo et al.l ([2002h . Fsoft has a similar but negative 
value and thus Frl is close to zero, independent of the value 
of rsoft - 

This illustrates the point made before: the interior of 
the Roche lobe is a complicated and variable dynamical sys- 
tem, and it is difficult to define a simple spherical region 
of radius rcut, that dynamically belongs to the planet and 
could be neglected in the torque calculation. Any fast migra- 
tion calculation would be very sensitive to the choice of rcut 
(fD'Angelo et al. 2005 ), and a wrong choice can give system- 
atic errors in the torque calculation. It is therefore better to 
include all of the torques, even though this implies dealing 
with the region very close to the planet. 

5.1.2 Mass accumulation in the planet's vicinity 

Up to this point we studied models with a constant planet 
mass, neglecting the gravitational interaction of the planet's 
envelope with the whole disc. However, the mass of gas ac- 
cumulating in the planet's vicinity in the case of constant 
hp = 0.4 is comparable to the initial planet mass and can 
influence the orbital evolution. In this section we present 
models where we vary the planet's mass. 

To test the dependence of the planet's orbital evolution 
on a varying planet mass we performed three simulations. In 
the first one the effective planet mass is constant during the 
whole simulation Mp = Mp = 0.001, similar to the models 
presented above. In the second model the planet mass was 
replaced with all of the mass within the smoothing length 
rsoft around the planet (Mp = Mp). In the last model we 
studied the planet growth through gas accretion, i.e. remov- 
ing gas from the planet's environment and adding its mass 
and momentum to the planet (see Sect I2.1.3|) . 

The results are presented in Fig. Illj^Curves 1 and 2 cor- 
respond to the constant mass and the Mp case respectively. 
Curve 3 presents the accreting model. The upper left and 
upper right panels show the orbital evolution and the non- 
dimensional migration rate Z (see below). The planet mass 
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Figure 9. The planet orbital evolution for different smoothing lengths of the planet's gravitational field for /ip = 0.4. The time evolution 
of the planet's semi-major axis a (left panel) and the mass of the gas inside a Hill sphere (right panel) are plotted. The planet mass 
Mp = 0.001. Curves 1 and 2 correspond to rgoft equal 0.0208 and 0.3-Rh respectively. 
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Figure 11. Results of the simulations for the varying planet's mass and /ip = 0.4. Curve 1 corresponds to the model with the constant 
planet's mass Mp = 0.001. Curve 2 shows the simulation with the effective mass of the planet augmented with the mass within the 
planet's smoothing length (Mp). Curve 3 presents the model with an accreting planet. Upper left and upper right panels show the 
orbital evolution and the non-dimensional migration rate Z respectively. The planet mass and the mass of the gas inside a Hill sphere 
are presented on the lower left and lower right panels. During the first 70 orbits (fast migration limit |Z| > 1) all systems evolve almost 
exactly the same way. Although in the fast migration regime the migration rate a is identical in simulations 1 and 2 and the differences 
for the accreting model are small, Z differs slightly, since df depends on the planet mass. The simulations start to differ after reaching 
the slow migration limit. 
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Figure 12. Results of the simulations of inward migrating planet with the varying planet's mass for hp = 0.3. Curve 1 corresponds to the 
model with the constant planet's mass Mp = 0.001. Curve 2 shows the simulation with the effective planet mass increased by the mass 
within the planet's smoothing length. Upper left and upper right panels show the orbital evolution and the non-dimensional migration 
rate Z for the inward migrating planets respectively. The planet mass and the mass of the gas inside a Hill sphere are presented on the 
lower left and lower right panels. 



and the mass of the gas inside a Hill sphere are presented 
on the lower left and lower right panels. 

The non-dimensional migration rate Z is defined as the 
ratio of the migration rate d and the so-called fast migration 
speed hi 



Z = 



(21) 



where af is given by the ratio of the half width of the 
horseshoe region Xs and t he libration time-scale Tub. In 
a Keplerian disc df eq uals (|Artvmowicz Peplinskil l2007l , 
IPapaloizou et al.ll2QQ7f ) 



af 



(22) 



Svra 

where Xs is estimated to be about 2.5 Hill sphere radii. We 
can divide the type HI migration into a fast, \Z\ > 1, and 
a slow, \Z\ < 1, migration regime. Type II- like migration 
phase corresponds to \Z\ <^ 1. The use of Z will be discussed 
more fully in Paper II. 

During the fast migration phase (lasting for about 70 
orbits) the planet's orbital evolution is almost independent 
of the planet's mass and the mass of the circumplanetary 
disc. All three models show a similar evolution even though 
the planet's effective mass can differ up to 40% and the 



mass of the circumplanetary disc can differ up to 100%. The 
amount of the gas flowing into the Hill sphere (for the model 
with accretion the sum of the Hill sphere content and the 
increase of the planet's mass) is similar too. This means 
that the inflow into the circumplanetary disc is a dynamical 
effect that mostly depends on the value of Z and the initial 
disc surface density, and only weakly depends on the gas 
dynamics deep in the Roche lobe. On the other hand, it 
should be very sensitive to the temperature profile at the 
boundary of the Roche lobe. Comparing the lower left and 
the lower right panels of Fig. [11] we can see that most of 
the mass within the Roche sphere is contained within the 
planet's smoothing length. 

The differences between the three cases show up in the 
slow migration phase, where the first and the second model 
lose gas from the Hill sphere. The simulation with constant 
Mp keeps losing mass from the circumstellar disc during that 
entire phase, whereas in the second simulation the amount 
of mass within the Hill sphere converges to 1.22Mp. In the 
accreting model the amount of gas available in the planet's 
vicinity and the pressure gradient are too small to cause 
mass loss from the circumplanetary disc. Instead the gas 
orbiting the planet is accreted, and in the stage of a gap 
creation (when the strong mass inflow into the circumplane- 
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tary disc stops) the amount of the mass inside the Hill sphere 
quickly decreases. At the end of the simulation the planet 
has reached a mass of 1.4Mp. 

The differences in the slow migration regime (\Z\ < 1) 
lead to different final positions in the three cases. The rea- 
son is that ttf depends on the planet mass. The case of a 
constant planet mass makes the transition to the slow phase 
latest, and thus achieves^the smallest orbit at a — 1.1. The 
difference between the Mp case and the accreting planet is 
due to the fact that the former loses mass from the Hill 
sphere. This influences the migration during the end of the 
slow migration phase (see the difference between curves 2 
and 3 between 70 and 90 orbits in the upper right panel of 
Fig. Ill|) and allows the accreting planet to migrate further 
(a = 1.2) than the planet in the model with the effective 
mass increased (a = 1.3). 

As we can see, the orbital evolution differs between the 
fast (|Z| > 1) and the slow (\Z\ < 1) migration limit. During 
the fast migration under constant hp, we can neglect details 
of the gas evolution deep in the circumplanetary disc, since 
the rate of the mass accumulation in the planet's proximity 
is relatively low and the migration is driven by the outer part 
of the Roche lobe. However, we should keep in mind that a 
variation of hp can modify the flow in the whole Roche lobe 
and thus can influence migration. The gas evolution in the 
circumplanetary disc becomes important for \Z\ ^ 1, when 
the planet slows down and starts to open a gap, and the rate 
of the mass accumulation in the planet's proximity increases. 

To test this last dependence we performed simulations 
with constant planet mass (case 1 above) and with the 
planet's effective mass increased by the mass within the 
planet's smoothing length (case 2 above) for hp equal 0.3, 
0.5 and 0.6. The results for hp = 0.3 are presented in Fig. [12] 
(curves 1 and 2 for Mp* = 0.001 and Mp* = Mp respec- 
tively). The upper left and upper right panels show the or- 
bital evolution and the non-dimensional migration rate Z, 
respectively. The planet mass and the mass of the gas in- 
side a Hill sphere Mh are presented on the lower left and 
lower right panels. We present here the results for hp = 0.3 
even though the planet's migration is slightly dependent on 
Tsoft for this value of the circumplanetary disc aspect ratio 
{Mh grows from 0.0015 up to 0.0022 for rsoft changing form 
O.SRu to 0.0208, and the migration rate a differs by less than 
5% between both models), since it illustrates the interesting 
case, were the sound speed Cs at the boundary of the Roche 
lobe given by E0S2 is smaller than Cs given by EOSl. The 
model with hp =0.4 is the limiting case, where both equa- 
tion EOSl and E0S2 give the same real disc aspect ratio 
outside the Hill sphere (see Fig.[T]). In the rest of the models 
E0S2 gives bigger Cs in the planet's vicinity than EOSl. 

The amount of mass accumulated in the planet's vicin- 
ity is very sensitive to the value of hp . UnUke the model with 
hp — 0.4, the model with /ip = 0.3 and Mp* — Mp allows 
the accumulation in the circumplanetary disc of an amount 
of mass comparable to the planet mass (about 50% of Mp 
after 10 orbits). In the second simulation (Mp* = Mp) this 
mass grows faster, since the Hill radius grows with Mp and 
the non-dimensional migration rate \Z\ decreases below 1, 
allowing a speed up of the accumulation (see Fig. [TT]) . Fi- 
nally Mp reaches 4 Jupiter masses and the planet migrates 
in the slow migration limit, stopping in the type II like mi- 



gration at about a = 2.2. As we can see th£ planet's migra- 
tion is sensitive to the rapid increase of Mp. The limiting 
case hp = 0.4 was already discussed and shows the differ- 
ent behaviour in the slow and the fast migration limit. In- 
creasing hp above 0.4 means significant decreasing Msoft and 
the planet's migration becomes independent of the choice of 
Mp*. This allows the planet to travel faster and further. A 
more detailed description will be presented in Paper II. 



5.2 Outward migration case 

Type HI migration has no predefined direction and can re- 
sult in both inward and outward directed migration. In the 
previous section we discussed the inward migration case. 
Most of the conclusions reached there are also valid for the 
outward directed migration, even though there are clear dif- 
ferences between the two types of simulations. 

The first important difference is that the amount of 
mass accumulated in the circumplanetary disc is much larger 
in the case of outward migration, up to several times the 
initial planet mass. This is caused by the combination of 
a different gas density at the planet's initial position, and 
the fact that the Roche lobe size (and the relative grid res- 
olution) grows in the case of of outward migration. In this 
section we focus on the dependence of the planet's migration 
on the value of the planet's gravitational softening and the 
treatment of Mp. 

5.2.1 Gravitational softening 

In Sect ion [5. 1.1 1 we presented the dependence of the planet's 
migration on rsoft for the inward migration case. We per- 
formed similar tests for the outward migrating planet. Fig- 
ure [13] shows the planet's orbital evolution and gas mass 
inside the Hill sphere for two simulations with rsoft equal 
0.0208 and O.SRu respectively. In both simulations we used 
hp = 0.4, and M^^Mp. 

In both cases the planet initially migrates outward but 
after about 35 orbits the direction of migration suddenly 
changes. We find that the planet's orbital evolution is inde- 
pendent of rsoft during the first 25 orbits. When the migra- 
tion speed starts to slow down, some differences between the 
two cases appear, and after the migration direction has re- 
versed a clear dependence on the smoothing length becomes 
apparent. Still, the migration rate in both simulations differs 
by less than than 20% (when comparing a at a given radius 
rather than at given time) and the mass of the gas inside 
the Hill sphere remains similar up to 60 orbits. The origin 
of the differences is an instability of the horseshoe region 
that arises for relatively high values of the migration rate. 
While the details of the instability are not entirely clear 
(streaming instability of gas with different vorticity or al- 
ternatively accelerating/decelerating migration of a planet 
might be involved), the net result is a decrease of the mass 
deficit (density contrast) of the librating disk region, which 
decreases the corotational torque. During the migration re- 
versal, the numerically computed migration rate depends 
strongly on rsoft • Therefore we limit our investigation in the 
case of Mp = Mp to the outward migration phase only, 
which is independent of rsoft- 

The dependence of the torques Fsoft, Tout and Frl on 
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Figure 13. Orbital evolution of the outward migrating planet for different smoothing lengths of the planet's gravitational potential. The 
time evolution of the planet's semi-major axis a (left panel) and the mass of the gas inside a Hill sphere (right panel) are plotted. 
Curves 1, 2 correspond to rgoft equal 0.0208 and O.SRu for hp = 0.4. The planet mass Mp = 0.001. 
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Figure 14. Results of the simulations of outward migrating planet with the varying planet's mass. Curve 1 corresponds to the model 
with the constant planet's mass Mp = 0.001. Curve 2 shows the simulation with the effective planet mass increased by the mass within 
the planet's smoothing length. Upper left and upper right panels show the orbital evolution and the non-dimensional migration rate Z 
for the inward migrating planets respectively. The planet mass and the^ass of the gas inside a Hill sphere are presented on the lower 
left and lower right panels. The migration is sensitive to the value of Mp. In the first simulation the phases of the outward and then 
inward directed rapid migration are visible, but the second models shows only the first phase. Instead the planet is locked in the relatively 
slow migration mode at around a = 1.7. However, during the first 20 orbits both systems evolve in the same way, even the the mass 
accumulated in the planet's proximity differ by a factor of three. 
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Tsoft is similar for both the inward and the outward migra- 
tion. 

5.2.2 Mass accumulation in the planet's vicinity 

In Section r5.1.2l we presented models for an inward migrating 
planet where we used different prescriptions for Mp. In the 
case of the outward migrating planet we performed similar 
simulations with rsoft = O.SRu, hp = 0.4 and /xd = 0.01. 
Here wejDnly consider the cases of Mp = Mp = 0.001 and 
Mp = Mp, i.e. we do not consider the accretion case. We 
present the resultsjn Fig. [TU Curves 1 and 2 correspond to 
the Mp and the Mp case respectively. The upper left and 
upper right panels show the orbital evolution and the non- 
dimensional migration rate Z, the planet mass and the mass 
of the gas inside the Hill sphere are presented on the lower 
left and lower right panels. 

The orbital evolution of the outward migrating planet is 
sensitive to the choice for Mp. In the Mp = Mp simulation 
the planet reaches a = 2.25 and rapidly changes its direction 
of migration (as above). For the Mp = Mp case the mass 
accumulated in the planet's proximity grows from an initial 
value of about 4.5 (Mp = 0.001) up to about 10.5 Jupiter 
masses. This behaviour is similar to the case of inward mi- 
grating planet with hp = 0.3. The increase of Msoft causes 
the planet to stop around a = 1.75 after 32 orbits, since 
the disc is not massive enough to support rapid migration 
of such a high mass object. 

The effects of neglecting the gravitational interaction 
between the massive planet's envelope and the whole disc 
are visible on the plot of the non-dimensional migration rate 
Z. In the Mp = Mp simulation the planet migrates in a fast 
migration limit with Z > 1 and reaches Z — 4 before it 
changes its direction of migration. For Mp = Mp Z stays 
below 1 and the planet migrates in the slow migration limit. 
However, the big difference in Z is caused mostly by the 
different planet's effective mass. During first 20 orbits the 
semi-major axis a differs by less than 7% and the migration 
rate d differs at most by 25%. During this initial period 
both systems evolve in the same way, even though the mass 
accumulated in the planet's proximity differs by a factor of 
three (about 2 and 6 Mq^ for the two models after 20 orbits). 
This is caused by the fact that for Z > 1 the migration 
rate depends only weakly on the planet's mass, and during 
the first 20 orbits Z is close to Z = 1 in the Mp = Mp 
simulation. ^ 

This test shows that it is important to use Mp when we 
want to study the slow migration limit of migration type HI 
and its stopping. However, just using a constant Mp gives 
a correct description of the fast migration limit. For hp = 
0.6 the amount of mass accumulated in the circumplanetary 
disc and the planet's migration become independent on the 
treatment of the effective mass. A more detailed description 
of the outward- directed type HI migration will be given in 
Paper HI. 



6 CONCLUSIONS 

We investigated the orbital evolution of a Jupiter-mass 
planet embedded in a disc and the dependence of its migra- 
tion on the adopted disc model, as well as the representation 



of the disk-planet coupling. The calculations were performed 
in two dimensions in an inert ial frame of reference, using a 
Cartesian coordinate system. We focused especially on the 
details of the flow in the planet's proximity. We employed an 
adaptive mesh refinement code with up to 5 levels of refine- 
ment which provides high resolution in the planet's vicinity. 
We conducted a careful analysis of the effects of various nu- 
merical parameters that influence the flow patterns inside 
the Roche lobe and the dependence on the grid resolution. 
Our goal was to find a disc model that allows numerical 
convergence to be reached and that minimises non-physical 
effects such as an artificial increase of the planet's inertia. 
We introduced a first-order correction of the gas accelera- 
tion in the circumplanetary disc due to the gas self-gravity 
and a modification of the local isothermal approximation, 
which allows us to take into account the gravitational field 
of the planet when calculating the local sound speed. 

For a simple disc model without these modifications 
to the temperature profile and the gas acceleration, we 
obtain ed results consistent with those of iD'Angelo et al.l 
(|2005l ). depending strongly on resolution and gravitational 
smoothing. Convergence could only be achieved if the 
torques from the interior of the Hill sphere Frl are ne- 
glected. However, the interior of the Roche lobe is important 
for driving rapid migration and cannot be neglected when 
numerically studying type HI migration. Since this region is 
complex and dynamically variable, it is impossible to define 
a unique radius around the planet that would separate the 
planet's sphere of influence from the disk, and thus could 
be neglected in the torque calculation. Therefore, one has to 
include all the torques from the Roche lobe, even though ar- 
tiflcial gravitational softening must be employed on spatial 
scale which is much smaller than the Roche lobe. Moreover, 
in the case of fast migration, the interior of the Roche lobe 
is not separated from the corotational flow and we cannot 
neglect the gas inflow into the circumplanetary disc. 

When including the torques from the entire Roche lobe, 
we achieve convergence when using the modiflcation of the 
local isothermal approximation and the correction for self- 
gravity. Both corrections are needed and we need to use a 
scale height of the circumplanetary disc hp ^ 0.4 to achieve 
this convergence. In this case the results are insensitive to 
the gravitational smoothing length rsoft- 

However, the (converged) solution will depend on the 
choice for hp. It determines the flow structure inside the 
Roche lobe and the amount of gas accumulated in the 
planet's vicinity. The latter influences the migration be- 
haviour as the accumulation of a large mass lowers the non- 
dimensional migration rate Z. As long as \Z\ > 1, the planet 
migrates rapidly at a rate that is only weakly dependent on 
the planet's mass and the conditions in the circumplanetary 
disc. If the mass accumulation forces |Z| to drop below 1, mi- 
gration will slow down and the planet may be locked in the 
disc in a mode resembling type II migration where, strictly 
speaking, corotational torques may still be dominant over 
Lindblad torques. 

Our modiflcations of the disc model compared with pre- 
vious investigations are physically justified, since the rapid 
accretion during migration will release potential energy and 
heat the gas in the planet's environment, while self- gravity 
will start to play a role if the planet collects an amount of 
mass of order its original mass. However, our implement a- 
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tion of these two effects are admittedly still very crude ap- 
proximations. Ideally, non-isothermal, self- gravitating mod- 
els should be used to self-consistently study type III migra- 
tion. For example, the assumption that the circumplanetary 
disk aspect ratio hp is constant during the whole simulation 
is probably not realistic. However, type III migration is only 
weakly dependent on the mass accumulation in the fast mi- 
gration Umit (|Z| > 1), so we expect variations of hp to play 
an important role only in the slow migration limit (|Z| < 1). 

In summary, the migration rates of the freely migrating 
planet in a massive disc (satisfying the condition for rapid 
migration) are sensitive to the adopted disc model, espe- 
cially to the treatment of the interior of the Roche lobe. 
The effects of heating close to the planet and self-gravity of 
the gas have to be taken into account to construct consistent 
models for type III migration. These models then show rapid 
migration independent of the numerical resolution used. In 
Papers II and III we will explore the physics of type III 
migration in more detail. 
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APPENDIX A: NUMERICAL DIFFUSIVITY OF 
CODE 

Our code participated in the Hydrocode Comparison f or a 
Disc-Planet System Project (|de Val-Borro et aP l2006l l . In 
this project the results of a wide range of hydrodynamics 
codes were compared for the case of Jupiter- and Neptune- 
mass planets on a constant orbit. In this comparison differ- 
ences between our Cartesian and the other (mostly) cylin- 
drical codes were found, including a version of FLASH using 
a cylindrical grid. 

The main difference between the Cartesian and cylin- 
drical FLAS H (denoted (FLASH-AP and FLASH-AG re- 
spectively in (|de Val-Borro et aI1l2006h l is a depletion of the 
inner disc. This leads to a deeper and broader gap, an asym- 
metry between the tadpole regions around the L4 and L5 
Lagrangian points and a difference in the differential Lind- 
blad torque calculation. The depletion is an effect of mass 
accumulation in the star's vicinity, which takes place in a 
poorly resolved part of the disc (compared to the cylindri- 
cal codes). High er reso lution simulations (not discussed in 
Ide Val-Borro et al.ll2006l ) show that this behaviour depends 
strongly on the grid resolution. 

To further investigate this feature and its relevance for 
the simulations presented in the paper, we performed a test 
employing the typical mesh size adopted in our investiga- 
tion. We used as = and the Jupiter- mass planet is kept 
on a constant orbit with semi-major axis a — 1. The disc 
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Figure Al. The normalised surface density profiles averaged az- 
imuthally over 27r for the Jupiter-mass planet. The star's neigh- 
bourhood is not p l otted to make the plot similar to Fig. 4 in 
Ide Val-Borro et al.l (l2006l V The curves 1, 2 and 3 correspond to 
60, 150 and 300 orbits. 

extends from —2.0 to 2.0 in both directions around the cen- 
tre of mass and the grid has 800 cells in both directions. 
The normalised surface density profile averaged azimuthally 
over 27T for 60, 150 and 300 orbits are presented in Fig. lAll 
0. As type III migration for Jupiter in our simulations takes 
about 60 orbits on average, curve 1 is the most relevant. 
At the higher resolution the normalised surface density pro- 
file is almost ident i cal to the results of the other codes in 
Ide Val-Borro eTHI (120061 ). The depletion of the inner disk 
only appears at later times (curves 2 and 3) . This means that 
the gas accumulation in the star's vicinity is a slow process 
and does not play an important role for type III migration 
of a Jupiter- mass planet, which is anyway dominated by the 
co-orbital torques. Clearly for a study of the slower Type II 
migration, dominated by the Lindblad torques, this effect is 
less than beneficial. 



^ The results for Hydrocode Comparison for a Disc-Planet Sys- 
tem Project for different times and resolution are presented at 
|http: / / www.astro.su.se/groups/p lanets / comparison / codes.html | 
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